eegmp_analysis/tools/[MEG]PLS/MEGPLS_PIPELINE_v2.02b/MEGPLS_COMPONENTS/MEGpipeline_RemoveICAcomps.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reject selected ICA components and reconstruct MEG data. %
% Last modified: July 25, 2014                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2013-2014, Michael J. Cheung
%
% This file is a part of the MEG & PLS Pipeline (MEGPLS). For more 
% details, see the documentation included with the software package.
%
% MEGPLS is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License version 2 as published by 
% the Free Software Foundation. This program is distributed in the hope 
% that it will be useful, but WITHOUT ANY WARRANTY; without even the 
% implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program. If not, you can download the license here: 
% <http://www.gnu.org/licenses/old-licenses/gpl-2.0>.


function MEGpipeline_RemoveICAcomps(InputParams)


% Make sure java paths added:
UseProgBar = CheckParforJavaPaths;


% Load InputParams:
name		= InputParams.name;
paths		= InputParams.paths;
RejectComps = InputParams.RejectComps;


% Set up errorlog and diary:
if exist('ErrorLog_ICA.txt', 'file')
    system('rm ErrorLog_ICA.txt');
end
if exist('Diary_ICA.txt', 'file')
    system('rm Diary_ICA.txt');
end

diary Diary_ICA.txt
ErrLog = fopen('ErrorLog_ICA.txt', 'a');



%==============================%
% CREATE ICA CLEANED DATASETS: %
%==============================%

for g = 1:length(name.GroupID)
	NumSubj = length(name.SubjID{g});
    
    if UseProgBar == 1
        ppm = ParforProgMon...
            (['REMOVING COMPONENTS & RECONSTRUCTING DATA:  ',name.GroupID{g},'.  '], NumSubj, 1, 700, 80);
    end
	
	parfor s = 1:length(name.SubjID{g})
		for c = 1:length(name.CondID)
			
			UnmixingInfo = LoadFTmat(paths.ICAunmixing{g}{s,c}, 'ICA');
			if isempty(UnmixingInfo)
				continue;
			end
			
			
			%--- If no components selected for rejection: ---%
			%------------------------------------------------%
			
			% Keep data as is:
			if isempty(RejectComps{g}{s,c})
                ICAcleanMEG = LoadFTmat(paths.InputPreprocMEG{g}{s,c}, 'ICA');
                if isempty(ICAcleanMEG)
                    continue;
                end
                
				ICAcleanMEG.ICAcompRemoved = [];
			end
			
			
			%--- If components selected for rejection: ---%
			%---------------------------------------------%
			if ~isempty(RejectComps{g}{s,c})
                
                % Check input dataset:
                CheckInput = CheckPipelineMat(paths.InputPreprocMEG{g}{s,c}, 'ICA');
                if CheckInput == 0
                    continue;
                end
				
				% Decompose original data using unmixing info:
				cfgICA			 = [];
				cfgICA			 = UnmixingInfo.cfgICA;  % Use same settings
				cfgICA.unmixing  = UnmixingInfo.unmixing;
				cfgICA.topolabel = UnmixingInfo.topolabel;
                cfgICA.inputfile = paths.InputPreprocMEG{g}{s,c};
				
				OrigComp = ft_componentanalysis(cfgICA);
                
                % Make sure OrigData has same channels as specified in cfgICA:
                % Due to ft_rejectcomponent keeping all channels in.
                cfgSelectChan           = [];
                cfgSelectChan.channel   = cfgICA.channel;
                cfgSelectChan.inputfile = paths.InputPreprocMEG{g}{s,c};
                
                OrigData = ft_selectdata_new(cfgSelectChan);
				
				% Reject components and reconstruct data:
				cfgReject			= [];
				cfgReject.component = RejectComps{g}{s,c};
				cfgReject.demean    = cfgICA.demean;
				
				ICAcleanMEG = ft_rejectcomponent(cfgReject, OrigComp, OrigData);
				ICAcleanMEG.ICAcompRemoved = RejectComps{g}{s,c};
                
                % Make sure AllEventInfo field is carried over:
                if ~isfield(ICAcleanMEG, 'AllEventInfo') && isfield(OrigData, 'AllEventInfo')
                    ICAcleanMEG.AllEventInfo = OrigData.AllEventInfo;
                end
                
                OrigComp = [];  % Free memory immediately
                OrigData = [];
			end
			
			
			% Save ICA cleaned data:
			CheckSavePath(paths.ICAcleanData{g}{s,c}, 'ICA');
			ParforSaveData(paths.ICAcleanData{g}{s,c}, ICAcleanMEG);
            
            ICAcleanMEG = [];  % Free memory
			
		end  % Cond
        
        if UseProgBar == 1 && mod(s, 1) == 0
            ppm.increment();  % move up progress bar
        end
        
	end  % Subj
    if UseProgBar == 1
        ppm.delete();
    end
	
end  % Group



%=========================%
% CHECK FOR OUTPUT FILES: %
%=========================%

for g = 1:length(name.GroupID)
    for s = 1:length(name.SubjID{g})
        for c = 1:length(name.CondID)
			
			if ~exist(paths.ICAcleanData{g}{s,c}, 'file')
				fprintf(ErrLog, ['ERROR: Missing ICA cleaned dataset for:'...
					'\n %s \n\n'], paths.ICAcleanData{g}{s,c});
			end
            
        end
    end
end



%=================%

if exist([pwd,'/ErrorLog_ICA.txt'], 'file')
    LogCheck = dir('ErrorLog_ICA.txt');
    if LogCheck.bytes ~= 0  % File not empty
        open('ErrorLog_ICA.txt');
    else
        delete('ErrorLog_ICA.txt');
    end
end

fclose(ErrLog);
diary off

end



function ParforSaveData(OutputPath, ICAcleanMEG)
	ICAcleanMEG
	save(OutputPath, 'ICAcleanMEG');
end